#define VERBOSE

using System;
using System.Collections.Generic;
using System.Text;
using Algorithms.Numerical;

namespace Algorithms
{
    public class Invariants
    {
        /**
         *  Columns of resultant matrix form p-invariants. 
         *  This matrix can be degenerated (width 0) in case there is no p-invariants.
         */
        public static IntMatrix CalulatePInvariants(IntMatrix incidentionMatrix)
        {
            return CalculateInvariants(incidentionMatrix.Transponed);
        }

        /**
         * Columns of resultant matrix form t-invariants. 
         *  This matrix can be degenerated (width 0) in case there is no t-invariants.
         */
        public static IntMatrix CalulateTInvariants(IntMatrix incidentionMatrix)
        {
            return CalculateInvariants(incidentionMatrix);
        }


        private static IntMatrix CalculateInvariants(IntMatrix im)
        {
            Matrix m = im.Matrix;
            Matrix sbase = Algebra.Instance.SolveUniform(m, true);
            Matrix view = MakeNonnegative(sbase);

            int toSkipCnt = 0;
            int[] toSkip = new int[sbase.Width];
            for (int x = 0; x < sbase.Width; x++)
            {
                bool has = false;
                for (int y = 0; y < sbase.Height; y++)
                    if (sbase[y, x] < 0)
                        has = true;
                if (has)
                    toSkip[toSkipCnt++] = x;                    
            }

            view = sbase.SelectionView(null, selectionInverse(toSkip, toSkipCnt, sbase.Width));

            IntMatrix r = new IntMatrix(view.Height, view.Width);
            for (int a = 0; a < view.Height; a++)
                for (int b = 0; b < view.Width; b++)
                    r[a, b] = (int)view[a, b];
            return r;
        }

        protected static Matrix MakeNonnegative(Matrix mtx)
        {
            GaussMatrix g = Algebra.Instance.GaussEliminateBelow(mtx.TransposedView, true);
            g = Algebra.Instance.GaussEliminateAbove(g, true);
            g = Algebra.Instance.MakeIndepColumnsUniform(g);
            Matrix m = g;

            #if VERBOSE
            Console.WriteLine("MakeNonnegative():After gauss:");
            Console.Write(m.ToString());
            #endif
            
            int x, y, i;
            int plusIdx = -1;
            int[] minusIdxs = new int[m.Height];
            int minusCnt = 0;
            // for each column
            for (x = 0; x < m.Width; x++)
            {
                minusCnt = 0;
                plusIdx = -1;
                // for each row in column
                for (y = 0; y < m.Height; y++)
                    if (m[y, x] > 0)
                        // remember one positive element
                        plusIdx = y;
                    else if (m[y, x] < 0)
                        // remember all negative elements
                        minusIdxs[minusCnt++] = y;

                if (minusCnt > 0)
                {
                    if (plusIdx >= 0)
                    {
                        for (i = 0; i < minusCnt; i++)
                        {
                            int minusIdx = minusIdxs[i];

                            int gcd = Gcd( (int)m[minusIdx,x], (int)m[plusIdx,x] );
                            int plusMul = Math.Abs((int)m[minusIdx, x] / gcd);
                            int minusMul = Math.Abs((int)m[plusIdx, x] / gcd);

                            Matrix aux = m.SelectionView(new int[] { plusIdx }, null).Copy();
                            aux.Mul(plusMul);
                            m.SelectionView(new int[] { minusIdx }, null).Mul(minusMul);
                            m.SelectionView(new int[] { minusIdx }, null).Add(aux);
                        }

#if VERBOSE
                        Console.WriteLine("MakeNonnegative(): After column {0} fix", x);
                        Console.Write(m.ToString());
#endif

                    }
                    else
                    {
                        // All rows with negative elements in this row cannot be fixed. 
                        m = m.SelectionView(selectionInverse(minusIdxs, minusCnt, m.Height), null);
#if VERBOSE
                        Console.WriteLine( "MakeNonnegative(): After row skip. Col {0}", x );
                        Console.Write( m.ToString());
#endif
                    }
                }
            }

            return mtx;
        }

        protected static int[] selectionInverse(int[] sel, int length, int max)
        {
            if (length > 0 && sel[length - 1] >= max)
                max = sel[length - 1] + 1;
            int []inv = new int[max - length];
            int last = -1;
            int k = 0;
            for (int i = 0; i < length; i++)
            {
                for (int j = last + 1; j < sel[i]; j++)
                    inv[k++] = j;
                last = sel[i];
            }
            for (int j = last + 1; j < max; j++)
                inv[k++] = j;
            return inv;                
        }

        protected static int[] selection(int[] sel, int length)
        {
            int[] s = new int[length];
            for (int j = 0; j < length; j++)
                s[j] = sel[j];
            return s;
        }


        private static int Gcd(int a, int b)
        {
            a = Math.Abs(a);
            b = Math.Abs(b);
            while (b != 0)
            {
                int t = b;
                b = a % b;
                a = t;
            }
            return a;
        }

    }
}
